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SIMIARY 


A detailed description is given of how the decoupling approximation known 
as the doubly asymptotic approximation (DAA) can be implemented with NASTRAN to 
solve shock problems for submerged structures. The general approach involves 
locating the nonsymmetric terms (which couple structural and fluid variables) 
on the right-hand side of the equations. Ihis approach results in coefficient 
matrices of acceptable bandwidth but degrades numerical stability, requiring 
a smaller time step size than would otherwise be used. It is also shown how 
the structure’s added (virtual) mass matrix, a necessary ingredient to DAA, can 
be calculated with NASTRAN. The version of NASTRAN used is NASA’s standard 
level 16 with one program modification, velocity-dependent nonlinear loads, for 
which the FORTRAN changes are listed. 


STATEMENT OF THE PROBLEM 


The general class of problems known as fluid-structure interaction 
problems includes the special case of determining the shock response of sub- 
merged structures. This is of particular interest to naval engineers concerned 
with the dynamic structural response of submarines (including the hull, 
appendages , and internal equipment) . 

Consider the idealized situation consisting of a ring-stiffened cylindrical 
shell with flat end caps which is deeply submerged in water, initially at rest, 
and subjected to the shock of a distant underwater explosion (fig. 1) . The 
general problem is to conpute the time -dependent elastic structural response of 
the cylinder. We will further sinplify the problem with the following 
assumptions : 

1. The shock wavefront in the vicinity of the cylinder is planar, 
a reasonable assumption whenever the source of the shock is located far away. 

2. The time history of the free-field incident pressure is a step 
function. This assumption can be made without loss in generality since the 
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structural response to an arbitrary time variation of the pressure history can 
be easily obtained by a convolution integral involving the step response. (See 
Appendix D.) 

3. The shock wave hits the cylinder from the side rather than from 
the end or from some more general oblique direction. This assumption results 
primarily in simplifying the data deck which describes the problem. Removing 
the restriction presents no conceptual difficulty. 


BACKGROUND 


In either the absence or the presence of the surrounding water, the 
structure can be modeled with finite elements in the usual way: with plate or 
shell elements for the unstiffened cylinder, and with beam elements for the 
ring stiffeners. 

The fluid is generally treated as an acoustic medium (e.g., see ref. 1): 
a centres sible, inviscid fluid which undergoes only small amplitude motion and 
whose pressure p satisfies the wave equation 

v2p = p/c2 (1) 

where dots indicate time differentiation and c is the speed of sound. At an 
interface between the fluid and a solid 


9n 



( 2 ) 


where n is the unit outward normal from the solid at the interface, p is the 
fluid mass density, and u is the outward normal component of displacement of 
the interface . " 


In principle, the fluid part of the problem can be handled by modeling a 
portion of the fluid with finite elements (refs. 2, 3). In reference 3, for 
example, the analogy was drawn between the scalar wave equation (1) and the 
elasticity equations so that a standard structural analysis computer program 
like NASTRAN (refs. 4, 5) can be used to solve problems involving the wave 
equation, Poisson's equation, or Laplace's equation. For finite fluid regions 
such an approach presents no significant problems. However, for structures 
submerged in infinite fluids, the analyst is faced with the additional problem 
of truncating the fluid and applying a radiation condition at the artificial 
boundap^ in order to absorb outgoing waves. Even if a reasonable radiation 
condition could be formulated, the cost of explicit fluid modeling could be 
prohibitive . 

An attractive alternative to such modeling is provided by approximations 
(refs. 6, '') which uncouple the structural response from the fluid response in 
the sense that the fluid pressure at the fluid -structure interface is deter- 
mined (approximately) from a knowledge of only the interface motion. Although 
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several such decoupling approximations have been formulated, one which is 
currently attracting attention is the doubly asymptotic approximation (DM) of 
Geers (ref. 6) . This paper shows how the DM scheme can be cast in a form 
which can be conveniently solved by NASTRAN. The version of NASTRAN used is 
NASA's standard level 16 with one program modification: velocity- dependent non- 
linear loads. The less efficient earlier version of NASTRAN, level 15, can 
also be used, although at greater conputer expense. The addition of velocity- 
dependent nonlinear load capability is convenient, but not crucial, for the 
NASTRAN inplementation of DM. 


THE DOUBLY ASYMPTOTIC APPROXIMATION (DM) 


A submerged structure subjected to an underwater shock wave experiences, 
at any given time t, a total dynamic pressure p which can be considered to 
comprise two conponents: an incident pressure p^ which would occur if no 
obstacle (the structure) were present, and a scattered pressure pg which is the 
difference between the total pressure and the incident pressure. Thus, 

P = Pi + Pg C3) 

The scattered pressure p^ is sometimes further decomposed as 

Ps ' Prs * Pr 

where p is the scattered pressure which would result if the structure were 
rigid ana stationary and p^., the radiated pressure, is the remainder. Of the 
three components of pressure, only p;^. depends on the structural motion, 
whereas both pi and p^s can be computed as if the structure were rigid and 
stationary. 

For the submerged ring- stiffened cylinder of interest (fig. 1), the plane 
wave incident pressure Pi is taken to be a step function (with wavefront moving 
to the left) : 


Pi(x,t) = PqH(x - Xq + ct) (5) 

where H is the dimensionless Heaviside unit step function (zero for negative 
argument and unity for positive argument) , Xq-ct is the location of the wave- 
front at time t, c is the speed of sound, and pq is a constant. 

The scattered pressure pg, which depends on the structural motion and 
hence cannot be precomputed as a function of time, is determined by the doubly 
asymptotic approximation (DM) (ref. 6) from 

vdiere Pg is the vector of unknown scattered pressures at the wet grid points of 
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the structure, is the (full) added mass matrix for the structure (see 
Appendix A) , p and c are the fluid mass density and sound speed, respectively, 
A is a diagonal area matrix converting grid point pressures to grid point 
forces, and Us is the vector of scattered wave particle accelerations normal to 
the structure's surface. The bar is used to distinguish this vector from the 
complete acceleration vector iig, which involves all structural degrees of 
freedom rather than just the normal components at wet points. Surface normals 
are taken as positive going into the fluid. 

Surface normal accelerations, like pressures, are decomposed into incident 
and scattered components; hence. 


U 5 = u - u. (7) 

where u is the vector of total normal accelerations at the wet grid points, and 
U4 is the vector of normal components (positive into the fluid) of incident 
fluid particle accelerations. 

Equation (6) was designated "doubly asymptotic" because it exhibits the 
correct asymptotic behavior at both the low and high frequency limits : at the 
low frequency limit (vdiich normally corresponds to late time behavior for 
transient situations) , the first tern of (6) is dominated by the second term, 
and (6) reduces to 

Fs = Ap^ . «a“s (8) 

in which the fluid loading is due to added (virtual) mass effects alone. At the 
high frequency limit (early time behavior) , the first term of (6) dominates the 
second term, and (6) reduces to 

Ps = pcu^ (9) 

which is the usual radiation damping relation. Equations (8) and (9) are 
referred to individually as the virtual mass and plane wave approximations, 
respectively. 

In general, the DAA, equation (6), yields better results than either of 
the special cases, equations (8) or (9). Huang (ref. 8) compared a DAA 
solution to an exact solution for a spherical shell and found that the DAA 
solution had slightly faster oscillations and stronger dairying. Nevertheless, 
the DAA provides a good compromise between cost and accuracy for underwater 
shock problems. 


DAA WITH NASTRAN 


The differential equation of motion for the ring-stiffened cylinder of 
interest (fig. 1) is 
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( 10 ) 


Mii + Ku= - Ap= - A(p^ + p^) 


where u is the vector o£ unknown displacements at the grid points, M and K are 
the structure's finite element mass and stiffness matrices, respectively, p is 
the vector (of dimension equal to the number of wet grid points) of fluid 
pressures at_the wet grid points, and A is an area matrix converting pressures 
to forces. A is not square (and hence not diagonal) because the vectors u and 
p are of different dimension. A contains non-zeros (equal to the area contri- 
butions) only at the intersections of rows corresponding to wet structural 
degrees of freedom with the columns of associated pressure variables. Thus, each 
row and column of A has at most one non- zero entry. A reduces to A if the zero 
rows are deleted and if the wet structural degrees of freedom are sequenced in 
the same order as the corresponding fluid pressure degrees of freedom. The 
area matrices defined here are "lumped" rather than "cons is tent". _ To switch to 
consistent loading, one need only change the area matrices A and A. 


In equation (10) , the total dynamic fluid pressure p is deconposed into 
incident and scattered pressures given by equations (5) and (6) , respectively. 
Since (6) is a differential equation, the complete problem involves solving (10) 
and (6) simultaneously, where the right-hand side of (6) is replaced by its 
equivalent from eqtiation (7) . 


The incident fluid particle normal acceleration vector is computed as 
follows: In general, the ratio of the pressure to the volume strain defines 

the bulk modulus k. Since k = p c^ for the acoustic fluid, we have, for a plane 
wave. 


P = 



P c" 


9x 


( 11 ) 


In particular, for the incident component. 


Pi = - P c 


9U . 

2 

3X 


( 12 ) 


From (5) and (12) , it follows that 


p. = - p c u . 

XI 


(13) 


where Uj^i is the x- component of incident , fluid particle velocity. The normal 
component of incident particle velocity u^ is 


Hence 


u- = u. cos 0 

1 IX 


p^ cos 0 = - p c u^ 


(14) 


(15) 
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Because as given by (5) is a step function, the u^ needed on the right- 
hand side of (o) is a Dirac delta. This problem can be avoided by defining a 
new variable q such that 

q = Pg (16) 

and time integrating equation (6) . Eqioations (10) and (6) then become 

Mii+Ku=-Ap. -Aq 

-1 \ . ( 17 ) 

Aq+ p c AM Aq = pcA u+ Ap. cos 6 
«■ 1 


where the second of equations (17) has also been multiplied by the area matrix 
A to symmetrize the coefficient matrices, and 

Au = u (18) 

In matrix format, these equations are 


o]itt* [“ Ji^l 



u 

j-Ap^ - 

_ 1 

q 

1 j pc A u + Ap^ cos 6 


which is the forni of the equations which NASTRAN uses. 


(19) 


It is interesting to observe that the new variable q defined by equation 
(16) is, in essence, the (scattered) velocity potential, since for an acoustic 
fluid the velocity potential 4) is related to the pressure p by (ref. 9) 

P = - P i> (20) 

Thus, as a consequence of trying to avoid the numerical problem of a Dirac 
delta, the fundamental unknown for the fluid is switched from the pressure to 
the velocity potential, thus returning to the established convention of fluid 
dynamicists . 


In equation (19), the unknowns u and q are defined using GRID cards. For 
the variables q, only one degree of freedom per point is retained. The usual 
finite element modeling of the structure yields M and K. The damping matrix is 
created by attaching dashpots (®AMPi) between each interface fluid point and 
groimd. The fluid matrix pcAMa^A can be assembled either by supplying it 
directly (on EMIG or IMI cards) or by letting NASTRAN compute it using an 
explicit finite element model of a portion of the fluid region. In Appendix A, 
it is shown that 
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p cAM'^A = cH 

3 > 


(21) 


where H is the fluid stiffiiess matrix condensed to the wet degrees of freedom. 
This condensation (using OMIT cards) is not necessary for the calculation but 
may result in a faster integration. The multiplicative constant c in the term 
cH, equation (21), can be automatically incorporated by setting the shear 
modulus on the MATl card equal to c rather than unity. 

The right-hand side of (19) consists of both time -dependent and velocity- 
dependent loads, which are supplied using TLOADl and NOLINl cards, respectively. 
(See Appendix C.) The input of the incident pressure is particularly simple 
since it is a step function. The input data can be further simplified by using 
the DELAY card to indicate that the incident wavefront (which is traveling at 
speed c) does not reach all points at the same time. 

In eqioation (19) , M and A are diagonal matrices , and K and H are positive 
definite and symmetric. K is also large and banded. H can be either large and 
banded, or small and full, depending on whether a static condensation (with 
CMIT's) is applied to it. 

The unloiowns in equation (19) are arranged so that the structural and 
fluid variables are uncoupled on the left-hand side, the only coupling 
occurring on the right. Thus the grid points should be sequenced to maintain 
the indicated partitioning and to give both K and H the smallest possible matrix 
wavefront (refs. 10, 11). 

The time step size needed to achieve numerical stability when the velocity- 
dependent terms are on the right was found to be about 1/10 of the transit time 
(the time required for a wave to travel one radius of the cylinder at speed c) . 


ALTERNATIVE DAA APPROACHES 


' Since the velocity-dependent loads in equation (19) are linear, they can 

be moved to the left-hand side and incorporated in the damping matrix. Symmetry 
I (but not positive definiteness) is then retained by dividing the second equation 
' in (19) by - pc. Since this formulation causes fluid-structural coupling on 
I the left, the unknowns have to be sequenced taking into account this new _ 
connectivity. Tliis approach is practical only if the fluid stiffness matrix H 
1 is not condensed but left large and banded, so that the overall system can be 

I made banded. Otherwise, the added mass matrix coupling causes non- zeros far off 

I the matrix diagonals. 

f 

The principal advantage in placing the velocity-dependent terms on the left 
, is numerical stability, so that a larger integration time step can be used. 

With those terras on the right, as in equation (19), the matrix bandwidth (and 
hence wavefront) is jailer, and the user has the option of condensing the fluid 
^ "stiffness" matrix cH into the smaller, but full, matrix pcAM ^A. 

a 
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Another possible way to formulate the DAA is to make use of the decon^jo- 
sition of scattered pressure pg into rigid body scattered and radiated compo- 
nents (eq. (4)). Since only the radiation pressure Px« depends on structural 
motion, the rigid body scattered pressure p^s can be preconputed and combined 
with the incident pressure pi* In that case, the DAA (eq. (6)) must supply 
only Pj., which satisfies 

F>r + P cM^^Ap^ = p cu^ (22) 

where the normal conponent of fluid particle acceleration at the fluid-solid 
interface is decomposed into 


u = Ui + Urs + (23) 

For rigid stationary structures, equation (23) simplifies to 

Ui - u,3 = 0 (24) 

at the interface, so that, in general, 

u = u^ (25) 

at the interface. Thus, (22) is equivalent to 

F>r + PcM'^Ap^ = pcu (26) 

The advantage of this general approach is that the rigid body scattered 
pressure p^g computed in advance to whatever accuracy one wants, so that 

the only approximation remaining involves the radiation pressure p^. The dis- 
advantage, however, is that the pre -calculation (a nontrivial one) has to be 
done at all. The decision of whether to use equation (6) or (26) also depends 
on the relative sizes of p^^-g and pj., since if p^ were small it would not make 
sense to compute it accurately. Unfortunately, the relative sizes are problem- 
dependent and hard to estimate. 


EXPLICIT FINITE ELEMENT FLUID MODELING 


The problem of computing the linear shock response of submerged structures 
can, in principle, be solved by explicit finite element modeling of a portion 
of the fluid volume. The purpose of this section is to formulate the problem 
sufficiently so that it can be solved by NASTRAN once the user has picked a 
suitable radiation condition to apply at the outer fluid boundary. 

The total dynamic fluid pressure satisfies the wave equation (1) in the 
field. This pressure can be decomposed into the sum of incident and scattered 
pressures, p^ and p^, as in equation (3). Since p- is defined to satisfy (1), 
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Pg must also satisfy the wave equation (1). At a fluid-solid interface, the 
boundary condition (2) becomes 




P u 


(27) 


where, for the finite cylinder of figure 1 subjected to a plane wave incident 
pressure, equation (5), we have 

3p. 3p. 



Thus, frcHn (27), 

1 # = - ^ Pi cos 6 - p U (30) 

The above Neumann boundary condition is equivalent to specifying a "load" on 
each interface pressure variable Pg equal to 

A( ^p^ cos 6 + p u ) (31) 

where A is the area associated with the interface point, so that the resulting 
finite element equations take the form (ref. 3) 

[%]!;!■[' 

where here q, defined as in equation (16), includes all fluid points, not just 
interface points. The area matrices also have to be redefined slightly to 
reflect the change in dimension of the vector q. The above foimu^ation is 
consistent with the definitions of fluid inertia Q and stiffness H given in 
Appendix A, which differ from the definitions of ref. 3 by a constant factor 
p c^. 

Equation (32) is complete except for a radiation condition on the pressure 
variable q. Once the user decides what radiation condition to use, it can be 
incorporated into the matrix equation (32) . 

It is interesting to observe the similarity between the explicit finite 
element formulation, equation (32), and that which arises from the doubly 
asymptotic approximation, equation (19) . The right-hand sides and the overall 
stiffiiess matrices are the same in both cases . In (32) , the overall mass 


_ ul -AP. -A, 
cHj (Q) Apj^cos6+ pcA u 
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matrix now includes the fluid inertia, which NASTRAN computes whenever the user 
supplies a non- zero mass density, in this case equal numerically to 1/c. The 
damping matrix in (19) also appears in (32) if a radiation condition involving 
dashpots is used, although in (32) the dashpots connect the outer boundary 
points, rather than the interface pressure points, to ground. 

Another approach using explicit fluid modeling was recently described by 
Newton and Atchison (ref. 12) , who elected to use the full fluid pressure p 
(rather than pg) as the fundamental pressure unknown. In that case, the time- 
dependent part of the right-hand side of the equations of motion is replaced by 
a non- zero initial condition on p and p throughout the fluid region. 

The main inpediments to solving the shock problem by these approaches are 
the potentially high cost of modeling a three-dimensional region of fluid and 
the difficulty in determining the radiation condition. For one -dimensional 
problems, the correct radiation condition merely involves attaching grounded 
dashpots to the outer fluid boundary (ref. 2) . However, for general three- 
dimensional situations, the mathematically exact radiation condition is a more 
conplicated relation (vhich cannot be modeled using only masses, springs, and 
dashpots) coupling all pressure variables at the outer boundary (ref. 13) . 

Since the implementation of such a condition is impractical, the analyst must 
resort to approximate radiation conditions which will not absorb 100% of out- 
going waves. It is for reasons like these that decoupling approximations such 
as DAA are being used. 


REFERENCES 

1. Geers, T.L.: "Transient Response Analysis of Submerged Structures," Finite 

Eleme nt Analysis of Transient Nonlinear Structural Behavior, AMD-Vol. 14, 
edited by T. Belytschko, J.R. Osias, and P.V. Marcal, The American 
Society of Mechanical Engineers, New York, 1975, pp. 59-84. 

2. Zienkiewicz, O.C., and Newton, R.E.: "Coupled Vibrations of a Structure 

Submerged in a Compressible Fluid," Proc. Int. Symp. on Finite Element 
Techniques , Stuttgart, 1969. 

3. Everstine, G.C., Schroeder, E.A. , and Marcus, M.S.: "The Dynamic Analysis 

of Submerged Structures," NASTRAN; Users' Experiences, NASA TM X-3278, 
Sept. 1975, pp. 419-429. 

4. "The NASTRAN Theoretical Manual," NASA SP-221(03), Washington, D.C., 

March 1976. 

5. "TTie NASTRAN User's Manual," NASA SP-222(03), Washington, D.C., March 1976. 

6. Geers, T.L.: "Residual Potential and Approximate Methods for Three- 

Dimensional Fluid-Structure Interactions Problems," J. Acoust. Soc. Amer., 
vol. 49, no. 5 (part 2), 1971, pp. 1505-1510. 


216 


7. Clark, A. V. , Jr.: "A Study of Fluid-Structure Interaction and Decot 5 )ling 

Approximations," Naval Research Laboratory Report 7590, Washington, D.C., 
December 1973. 

8. Huang, H. : "A Qualitative Appraisal of the Doubly Asynptotic Approxima- 

, tion for Transient Analysis of Submerged Structures by Weak Shock Waves," 
NRL Memorandum Report 3135, Naval Research Laboratory, Washington, D.C., 
Sept. 1975. 

9. Newton, R.E. : "Finite Element Analysis of Two-Dimensional Added Mass and 

Danping," chap. 11 in Finite Elements in Fluids , vol. 1, ed. by 

R.H. Gallagher, J.T. Oden, C. Taylor, and O.C. Zienkiewicz, John Wiley and 

Sons, Ltd., London, 1975, pp. 219-232. 

10. Everstine, G.C.: "The BANDIT Conputer Program for the Reduction of Matrix 

Bandwidth for NASTRAN." NSRDC Report 3827, Naval Ship Research and 
Development Center, Bethesda, Md. , March 1972. 

11 . Everstine , G . C . : ' 'Recent Inprovements to BANDIT , ' ' NASTRAN: Users * 

Experiences , NASA TM X-3278, September 1975, pp. 511-521. 

12. Newton, R.E., and Atchison, D.L.: "Response of a Ring-Stiffened Cylinder 

to an Acoustic Blast Wave," Second International Symposium on Finite 
Element Methods in Flow Problems , S. Margherita Ligure, Italy, 

June 1976, pp. 701-713. 

13. Zarda, P.R. : "A Finite Element -Analytical Method for Modeling a Structure 

in an Infinite Fluid," NASTRAN: Users' Experiences , NASA IM X-3428 , 
October 1976. 

14. Protter, M.H. , and Weinberger, H.F. : Max unum Principles in Differential 

Equations , Prentice -Hall, Inc., Englewood Cliffs, N.J., 1967. 

15. Khabbaz, G.R. : "Dynamic Behavior of Liquids in Elastic Tanks," AIAA 

Journal , vol. 9, no. 10, October 1971, pp. 1985-1990. 

16. Gallagher, R.H. : Finite Element Analysis Fundamentals , Prentice-Hall, Inc., 

Englewood Cliffs, N.J., 1975. 

17. Streeter, V.L.: Fluid Dynamics, McGraw-Hill Book Company, Inc., New York, 

1948. 

18. Przemieniecki, J.S.: Theory of Matrix Structural Analysis , McGraw-Hill 

Book Company, Inc., New York, 1968. 

19. Glockner, P.G.: "Symmetry in Structural Mechanics," J . Struct . Div . , Proc. 

ASME, vol. 99, no. STl, Jan. 1973, pp. 71-89. 


217 


20. Joseph, J.A. , editor: "MSC/NASTRAN Application Manual for CDC 6000 

Series," MSR-32, The MacNeal-Schwendler Corporation, Los Angeles, 

Calif., 1975; also, pub. no. 86616700, Control Data Corporation, Data 
Services Publications, Minneapolis, Minn., 1975. 

21. Hildebrand, F.B.: Advanced Calculus for Applications. Prentice -Hall. Inc.. 

Englewood Cliffs, N.J., 1963. 


APPENDIX A - ADDED MASS MATRICES 


Consider an elastic structure submerged in a finite acoustic fluid, whose 
pressure p satisfies the wave equation 


V^p = p/c2 (Al) 

where c is the speed of sound in the fluid. If both structure and fluid are 
modeled with finite elements, the resulting matrix equations take the general 
form (refs. 2, 3) 


'm o' 

ju| . 

K A' 

-T 

L-pA qJ 

y 

_0 H_ 




(A2) 


where M and K are the usual structural mass and stiffness matrices, Q and H 
are the inertia and "stiffness" matrices for the fluid, A is the area matrix 
converting pressure to force at the fluid -structure interface nodes, and p is 
the fluid's mass density. 

In equation (A2) , H can be assembled from standard 3-D elasticity finite 
elements (ref. 3) if only the x-component of displacement at each point is 
retained (to represent the scalar quantity p) and Hooke's law is specified as 


(A3) 


In terms of the usual engineering constants, equation (A3) is equivalent 
(numerically) to choosing the shear modulus G and Young's modulus E as 


XX 


yy : 


zz 




yz 


xz ; 
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1 -1 -1 

-1 1 -1 

-1 -1 1 


XX 
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G = 1 


\ 

\ 

\ 


\ 

E = aG, a>>l 


whose a must be large enough so that 
from a. On most computers, a = 
three dimensions. In two dimensions 
are 


a+1 is indistinguishable (numerically) 
suffices. Equation (A4) applies only in 
(plane stress) , the corresponding constants 


G = 1 

E = BG, B<<1 


(A5) 


where B should not be so small that 1+B is indistinguishable (numerically) from 
unity. On most computers, B = 10*** suffices. 

Equation (AS) also applies to axisymmetric problems formulated in cylin- 
drical coordinates with axisymmetric elements such as NASTRAN's CTRAPRG. How- 
ever in this case only the z -component of displacement can be used to represent 
pressure, in contrast to Cartesian coordinates in which any of the three trans- 
lation components can be used. 

In equation (A2) , Q can be assembled from standard elasticity finite 
elements (ref. 3) if the mass density assigned to the material is numerically 
equal to l/c^. 

For an incompressible fluid, c->«> (or the frequency oj^O) and the wave 
equation (Al) reduces to Laplace's equation 

v2p = 0 (A6) 


j^so, Q = 0, so that p can be eliminated from (A2) to yield 

(M + p A H ^A^ ) u + K u = f 


(A7) 


thus defining the added mass matrix 

— ^ — -1 "T 

M _ p AH A 

3 . “ 


(A8) 


for the submerged structure. 

We observe that the area matrix A is ^non- square since the vectors u and p 
are of different dimension. In addition, A is such that each row or column has 
but one non- zero entry. This entry is the area assigned to a particular node 
and located at the row corresponding to the nodal outv;ard normal displacement 
and column corresponding to the associated pressure variable. Thu^, since A 
involves only the interface variables, the fluid stiffness matrix H in equation 
(A8) can be reduced by static condensation (Guyan reduction) prior to performing 
the matrix product in (A8). If the condensed stiffness matrix is denoted H, 
then the corresponding added mass matrix is 

II = pAH'^A (A9) 

a 
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where A is the diagonal area matrix and involves only wet degrees of 

freedom rather than all structural degrees of freedom as in M . 

d. 

A physical interpretation of a particular ij entry in or is that it 
is the normal fluid force induced at point i (on the fluid-structure interface) 
due to a unit normal acceleration at interface point j, with all other points 
held fixed. Thus, it is clear that is a fully populated matrix. Since the 
interface boundary condition is 


9n 



(AlO) 


the specification of acceleration at the interface is equivalent to a Neumann 
boundary condition. Hence, the calculation of the added mass matrix is 
mathematically equivalent to solving Laplace's equation (A6) in the fluid region 
with Neumann boundary conditions. For uniqueness, p must be specified somewhere 

This Neumann problem is also equivalent to a steady-state heat conduction 
problem in which one seeks the temperatures at all mesh points (on the fluid- 
solid interface) due to a unit heat source at each such point in turn. The 
matrix H in equation (A8) is exactly the heat conduction "stiffness" matrix 
computed by finite element heat transfer conputer programs if the thermal 
conductivity is specified as unity. 

Since, in heat conduction problems, the extreme temperatures must occur 
on the boundary, and uneven ten^erature distributions can be maintained only by 
supplying heat at the warmest point on the boundary and removing heat from the 
coolest point (ref. 14) , it follows that the individual elements of the added 
mass matrix are always positive. 

Thus far, this discussion of added mass matrices has assumed the fluid 
region to be finite. Of more interest in naval applications is the infinite 
region. In this case one can define and model a finite region of fluid whose 
outer boundary (where p=0) is "sufficiently far" from the structure. The major 
problem facing the analyst is deciding where to locate this outer boundary. For 
a given problem, one approach to insure that the outer boundary is distant 
enough is to compute the added mass matrix (condensed to include only wet 
degrees of freedom) with two different locations of the outer boundary and look 
for convergence of the dominant terms in the matrix. 

The calculation of added mass matrices for structures submerged in 
j^^finite fluids would be more appealing if it did not involve the explicit 
modeling of a portion of the fluid. Since the problem to be solved is a 
Neumann problem in the infinite region surrounding the structure, it can also 
be formulated in terms of simple sources distributed over the fluid-solid 
interface (ref. 15). For economy, the source density distribution is usually 
assumed constant over each surface element. Consequently all matrices 
(including the added mass matrix) refer to element centroids rather than to the 
finite element grid points. One possible approach for transforming an added 
mass matrix from element to grid point values is as follows: For simplicity 

assume a rectangular mesh of surface elements (fig. 2) , wliere a typical element 
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(number 4"^ is shown connecting grid points 5, 6, 15, and 16. The sin^lest 
relationship between the central and nodal displacements is the arithmetic 
average 

“4 " * “6 * “15 * “ 16 ^/“ 

A more complicated relation taking into account the element shape function can 
also be writteA. The complete transformation involving all elements is of the 
form 


A = ru (A12) 

where u is the displacement vector for the element centroids, u is the dis- 
placement vector for the grid points, and r is the transformation matrix. The 
added mass matrix can then be computed from 

M = M r (A13) 

which is the usual transformation relationship for finite element matrices 
(ref. 16) . In equation (A13) , is the added mass matrix referred to 
centroidal coordinates. The transformation (A13) may result in converting M^, 
which is non- singular, into a singular matrix M^. 


Virtual Mass 


It is of interest to relate the added mass matrix (as used here) to the 
added mass (virtual mass) defined by hydrodynamicists (e.g., ref. 17). Virtual 
mass is a scalar quantity defined and computed for rigid structures oscillating 
in a specific rigid body motion, e.g., heave of a ship hull form. Since the 
added mass matrix is general enough to allow arbitrary elastic structural 
motion, virtual mass is merely a special case. 

Recall that a physical interpretation of the added mass matrix is that a 
particular ij entry is the normal fluid force at point i due to a unit normal 
acceleration at j, with all other points held fixed. In computing virtual mass, 
the acceleration at all points is specified, and the component of force in a 
particular direction is desired. For example, if (j) is a vector describing the 
amplitude of the desired rigid body motion, the virtual mass m in the same 
direction is 


m = M (f> (A14) 

3 . 

where Mg^ is the added mass matrix and each component of <t> is equal to the 
cosine of the angle between the surface normal at a point and the direction of 
motion (assuming unit amplitude motion) . Equation (A14) is identical in form 
to the definition of generalized mass for vibration mode shapes (ref. 18) . 

In general, there exist six rigid body modes (three translations and three 
rotations), each of which induces six components of force. Thus (A14) can be 
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generalized to define a 6x6 rigid body virtual mass matrix m whose ij entry is 



M 


(A15) 


where is the vector specifying the rigid body acceleration and is the 
vector describing the direction of the generalized force induced. 


APPENDIX B - USE OF STRUCTURAL SYM^TRY 


In general, it is economically advantageous to exploit as much structural 
symmetry as possible when performing a structural analysis. This exploitation 
is possible whenever the structure, in the absence of loads, possesses geometri- 
cal and structural symmetry (ref. 19) . Since time -dependent nonsymmetric 
loads can always be decomposed into the sum of symmetric and anti-symmetric 
parts, the overall problem can be decomposed in the same way. The purpose of 
this appendix is to summarize how this decomposition works for the class of wave 
problems arising in computing submerged shock response. 

Consider the cylinder cross-section shown in figure 3 with a typical point 
number 1 and its image point number 2. The applied loads at the two points can 
be arbitrary functions of time. 

The cylinder possesses numerous planes of symmetry, including the xz- and 
yz-planes. Thus, only one-fourth of the circumference has to be modeled. (In 
this particular case, the structure is axisymmetric and thus can be modeled 
using axisymmetric elements with nonsymmetric loading.) The indicated loading 
is symmetric with respect to the xz -plane and nonsymmetric with respect to yz. 
Since the problem is linear, the loading can be decomposed as shovmi in figure 3. 
with 


F^(t) = (F^(t) + F2(t))/2 
FaW = (Fi(t) - F2(t))/2 


(Bl) 


where the decomposition results in one problem which is symmetric with respect 
to the yz-plane and another problem which is anti -symmetric with respect to yz. 
For each component part, it suffices to model but one quadrant (fig. 3) and 
apply the appropriate boundary conditions (either symmetric or anti-symmetric) 
for all points in the symmetiy planes. 

For structural grid points (whose fundamental unknown is displacement) 
lying in a plane of s>Tmnetry , the boundary conditions are that the points can 
suffer no translation out of the plane of symmetry and no rotation about in- 
plane lines. The anti -symmetry boundary conditions are that the complementary 
degrees of freedom are constrained. For example, in figure 3, all points 
lying in the yz-plane must satisfy 


222 


0 for symmetry 




0=0 for anti -symmetry 


(B2) 


where e denotes rotations. 

For fluid grid points (vhose fundamental unknown is pressure) lying in a 
plane of symmetry or anti -symmetry, the boundary conditions are 


= 0 for symmetry 
p = 0 for anti -symmetry 


(B3) 


In finite element analysis, the above symmetry condition on p is a natural 
boundary condition and thus automatically satisfied if the unknown p is left 
free. 


These conditions on fluid pressure grid points are applicable for any 
pressure points lying in a plane of symmetry or anti -symmetry, including- those 
occurring in an explicit modeling of the fluid volume for the purpose of 
computing added mass matrices; i.e., the added mass matrix also has to exhibit 
the proper symmetry. 


APPENDIX C - VELOCITY-DEPENDENT NONLINEAR LOADS 


The finite element formulation derived to implement the doubly asymptotic 
approximation (DAA) with NASTRAN involves loads which, at each time step, depend 
explicitly on the current structural motion rather than on time. The standard 
versions of NASTRAN (levels 15 and 16) currently allow displacement -dependent 
loads but provide no convenient way to specify loads which depend on velocity 
or acceleration. (The implementation of such loads with a combination of 
transfer function (TF) and nonlinear load (NOLINi) cards is not only inconven- 
ient but also results in nonsymmetric matrices.) 

NASTRAN can be easily modified to allow the user to apply velocity- and 
acceleration-dependent loads using the NOLINi cards now used only for dis- 
placement-dependent loads. This appendix summarizes the FORTRAN changes to 
NASTRAN (level 16) needed to implement such loads. 

The approach taken is compatible with that used in MSC/NASTRAN (ref. 20), 
in which the liser indicates velocity dependence by adding 10 to the displace- 
ment component number CJ or CK on the NOLINi card. This modification is 
extended here to allow acceleration dependence, which is indicated by adding 20 
to CJ or CK. Acceleration dependence as implemented here, however, is not fully 
general, since it does not allow a change in the time step size. The velocity 
dependence is fully general. 
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(Cl) 


The finite difference formulas lised to compute velocity and 
acceleration at the nth time step are 

“n ‘ ■ “n-l)/*' 

“n ' '“n ■ 2“n-l (^2) 

where u^^ is the displacement vector at the nth time step and At is the time step 
size. 

The listing of the FORTRAN changes to NASTRAN (level 16) appears in 
figure 4, in which the format of CDC's UPDATE utility is used. These modifica- 
tions were adapted from similar changes made to NASTRAN 's level 15 by Messrs. 
James M. McKee and N^^les M. Hurwitz of the David W. Taylor Naval Ship Research 
and Development Center. 


APPENDIX D - RESPONSE TO ARBITRARY TIME -DEPENDENT 
LOADING BY CONVOLUTION 


When the linear shock response of large conplex structures is to be 
conputed with NASTRAN, it is often preferable to conpute first the response to 
a step function, because (1) input data preparation for NASTRAN is simplified 
considerably, and (2) the response for any arbitrary time-dependent loading can 
be easily conputed later by a convolution (superposition) integral (e.g., 
ref. 21), the formulas for which are summarized here. 

Consider the general equation 


Lw(x,t) = f(t) 


(Dl) 


where L is a linear differential operator, w is some response variable (e.g., 
displacement, velocity, stress, etc.), and the forcing function f is considered 
here to represent the incident free- field pressure which arises in underwater 
shock problems. 


If w (x,t) is the response to a unit step function, then 


w 


(x,t) = f(0)w (x,t) + / f ( t)w (x,t-x)dT 


0 


f(t)w (x,0) + / f(x)w' (x,t-T)dx 
0 ^ 


(D2) 


where we define 


(D3) 
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Either of the quadrature formulas (D2) can be used to confute the response to 
an arbittary forcing function f (t) . Since the two relations give different 
results numerically, the convolution can be confuted both ways and averaged. 
For our work, a short conputer program was written to conpute w, given tabu- 
lated values of f and w for non-uniform spacing of the abscissas. 


225 




Figure 1. - Ring -Stiffened Cylindrical Shell 
Subjected to Plane Wave Shock 





1 



» 

5 

• 

• 

6 



• 

•4 

‘ 



15 

f 

• 

16 

• 




1 1 

1 

1 


• • « 


Figure 2. - Rectangular Mesh of Surface Elements 
for Added Mass Matrix Transformation 





Figure 3. - The Superposition o£ Symmetric and Anti -Symmetric 
Solutions for Time -Dependent Problems with 
Nonsymmetric Loading 
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*IDE*<T E>»42976 

•COHPILE 0 P 04 ,IFSlP,TR 01 CtTRan 
•OetETS 9P04.70,OP04,78 
0 

C VELOCITY- OR ACCELEPATXON-OSPtNOENT NL LOftO 
C 

C LOOK FOR COHPOMENTS .6T. 6 FOR IMOEPENOENT OOF ANO PACK WITH 
C NOLIN CARO TYPE 

C VELOCITY (OR ACCELERATION) IS INOtCATEO BY AOOIMG 10 (OR 20) TO 
C CONPONENT 

C 

IUOOT«0 

IF(eUF(6).LE.6) GO TO 1341 
IU0OT*(BUF(6f/10) * 10 
BUF(6)*BUF(G)-IU00T 
1341 IF(I1«NE*2) GO TO 134$ 

1F<9UF(6)»LE.6) GO TO 134$ 

XG£Ts(auF<a)/iin « lo 

eUF(8)s8UF(6)-IGET 
IUOOT*IUOOT410*IGET 
•DELETE 0PO4.87 

SUF(2) » II * lUOQT 
•DELETE 0PO4,163»OP04.164 

IGETs 3UF(1) - (BUF(1)/10I • 18 
IF(IG£T*NE*2«AN0.KK.E0.8I GO TO 1396 
•DELETE IFS1P.243 

TF( (H(6)/10) »6T.2I GO TO 8 
IF( (H(6I- (ii(6)/10l • 18) .GT.6) GO TO 8 
•DELETE IFS1P.250 

IF( <H(8)/10) *67.2) GO TO 8 
XF( (N(8)- (H(8)/10) • 10) .GT.6) GO TO 8 
•INSERT TR01C.28 

2 »IU3 
•INSERT TR010.13 

3 •1U1,0ELTAT,IU3 

C IN NL LOAD CALCULATION, lU IS LATEST U, XUl IS ONE 9ACK, AND 
C 1U3 IS TNO 3ACK 

•DELETE TR010.125 

IGETsIKL) -aZ(L)/10) • 13 
IFdGET.NE.l) GO TO 110 
•INSERT TROlO.iSi 
C 

C VELOCITY- ANO ACCELEPAT10N-0EP&N0ENT NL LOADS 

C 

C STRIP VELOCITY OR ACCEL FLAGS FRCH NOLIN CARO TYPE AMO COMPUTE 
C VEL OR ACCEL FOR FLAGGED COMPONENTS .... 

C OOOT « (O(N)-U(N-D) /OELTAT 

C UOOTOOT* (U(N) - 2 • U(N-l) ♦ U (N-2) ) /3ELTAT*^2 

C 

C 

C NOTE - - . ACCELERATION-OEPENOENT NL LOADS DO NOT MORK IF A 
C CHANGE IN TIME STEP SIZE OCCURS. 

C 

M«1.0/0ELTA T 

•OcLcTE TR010.185,TR01D.196 
MM«IU1*IZ(I^2) 

MMM«IU3«IZ(I*2) 

JFLG«IZ(I) /lot 

IFLG*(IZ(I)-JFLG*100)/10 

L*IZ(I)-JFLG^100-IFLG^10 

IFLG»IFLG+l 

jFLGsjFLG^l 

GO TO (192,193,194), IFLG 
C OISPLACEHENT-OEPENOENT 

192 XaZCll 

GO TO 196 

C VELOCITV-OEPENOENT 

193 X»(2(M)-Z(MM) )^M 
GO TO 196 

C ACCELERATION-OEPENOEMT 

194 XsA. 

IF(ICOUNT.GT.l) 

-X*(2(M)-2.*2(HM) »Z(MHN))^HH 
196 CONTINUE 
•DELETE TR010.197 

HN«1U1*IZ(I«4) 

MHM*IU3»IZ(It4) 

GO TO (212,213,214), JFLG 
C OISPLACEMENT-OEPENOENT 

212 FX*X*Z(N) 

GO TO 248 

C VELOCITY-OEPENOENT 

213 FXaX*<Z(M)-2(HM) )*H 
GO TO 248 

C ACCELERATION-OEPENOENT 

214 FX«0. 

IF(ILOOP.GT.l.OR.ICOUNT.GT.l) 

-FX*X*(Z(M)-2.*Z(MN)»ZIMNM) |*MH 


Figure 4. - Listing of FORTRAN Changes to NASTRAN Level 16 
for Velocity- and Acceleration-Dependent Nonlinear 

Loads 
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